Although this story looks as if it should be a perfect fit for an SEM, there are a few roadblocks: - percent cover and SGD are highly correlated/covaried. Need to figure out how to deal with that in an SEM because SEM does not accept separate models with the same response variable, so SGD and percent cover would need to be together, but then results are not accurate depiction of the situation.
The zeros are removed from Hyp 3 because only a few indivdual recruits on the plates and they were skewing the data
Considering Varari and Cabral in separate models; however, maybe combine all into one model and use Location as an interaction term
NOTE: needed to download and install an old version of “heavy” package because it was removed from the CRAN repository in 2022
## [1] 1
## [1] 146
## [1] 21.75
## [1] 0
## [1] 1
## [1] 0.5144263
## [1] 2
## [1] 131
## [1] 33.8
## [1] 0
## [1] 1
## [1] 0.3693079
Using a generalized linear model, negative binomial
Results: There is a site difference for the effect of silicate, which makes sense because the two sites are so different in their dynamics (silicate concentration, flow, etc). However, there is no effect of site for percent cover. Percent cover significantly affects the number of total recruits that settle on the tiles; however silicate only affects the total number of settlers at Varari (linearly and non-linearly), but not at Cabral. According to the AIC, the silicate model is a lower rank, so percent cover is a better/stronger driver of total settlement, but silicate is still important. Looking at the standardized effect sizes, the power of the relationship is strong with ….
Why separate it into 2 models? Because too many variables at play, and silicate & percent cover are too highly correlated
#######################
## decision is to keep both sites together and use "Location" as an interaction term in the model
#######################
############
#### Hypothesis 1a - percent cover
PocCover_negbinom <- glm.nb((sum_total+1) ~ log(percent_cover+1)*Location, data=SettlementDF_Full)
summary(PocCover_negbinom) ## RESULT: log percent cover is significant, no location interaction or just location effect
##
## Call:
## glm.nb(formula = (sum_total + 1) ~ log(percent_cover + 1) * Location,
## data = SettlementDF_Full, init.theta = 1.109358935, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.1524 0.3984 5.402 6.58e-08 ***
## log(percent_cover + 1) 0.7256 0.2827 2.566 0.0103 *
## LocationVarari 0.8417 0.5169 1.628 0.1035
## log(percent_cover + 1):LocationVarari -0.3824 0.3381 -1.131 0.2581
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.1094) family taken to be 1)
##
## Null deviance: 56.158 on 39 degrees of freedom
## Residual deviance: 44.170 on 36 degrees of freedom
## AIC: 349.05
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.109
## Std. Err.: 0.236
##
## 2 x log-likelihood: -339.049
#r.squaredGLMM(PocCover_negbinom)
## AIC = 345.31, p-value = 0.01
## try without logging percent cover
PocCover_negbinom2 <- glm.nb(sum_total ~ percent_cover*Location, data=SettlementDF_Full)
summary(PocCover_negbinom2) ## RESULT: slightly less significant, but same results
##
## Call:
## glm.nb(formula = sum_total ~ percent_cover * Location, data = SettlementDF_Full,
## init.theta = 0.9627769765, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.39356 0.34385 6.961 3.38e-12 ***
## percent_cover 0.18364 0.07744 2.371 0.0177 *
## LocationVarari 0.57558 0.45594 1.262 0.2068
## percent_cover:LocationVarari -0.11861 0.08345 -1.421 0.1552
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.9628) family taken to be 1)
##
## Null deviance: 56.298 on 39 degrees of freedom
## Residual deviance: 44.715 on 36 degrees of freedom
## AIC: 345.31
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.963
## Std. Err.: 0.204
##
## 2 x log-likelihood: -335.309
## AIC = 345.31, p-value = 0.01
check_model(PocCover_negbinom2)
############
#### Hypothesis 1b - SGD
SGD_negbinom <- glm.nb(sum_total ~ (silicate_avg + I(silicate_avg^2))*Location, data=SettlementDF_Full)
summary(SGD_negbinom) ## RESULT: silicate at Varari is significant both linear and nonlinear
##
## Call:
## glm.nb(formula = sum_total ~ (silicate_avg + I(silicate_avg^2)) *
## Location, data = SettlementDF_Full, init.theta = 0.9122803951,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.72899 2.37335 0.729 0.4663
## silicate_avg 0.62979 1.51219 0.416 0.6771
## I(silicate_avg^2) -0.06483 0.21840 -0.297 0.7666
## LocationVarari -5.10510 3.85202 -1.325 0.1851
## silicate_avg:LocationVarari 7.17745 3.63825 1.973 0.0485 *
## I(silicate_avg^2):LocationVarari -1.94437 0.85219 -2.282 0.0225 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.9123) family taken to be 1)
##
## Null deviance: 53.537 on 39 degrees of freedom
## Residual deviance: 45.038 on 34 degrees of freedom
## AIC: 351.87
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.912
## Std. Err.: 0.191
##
## 2 x log-likelihood: -337.872
## AIC = 351.39
check_model(SGD_negbinom) ## obviously very collinear bc linear and nonlinear
#######################
## plot the negative binomial
#######################
#######################
## standardized effect sizes
#######################
# Note: won't work on dependent var bc changes count data to continuous which then doesn't use glm.nb()
## only calculate sef for independent vars of SGD, SGD^2, and percent cover
## and take sef before the model, create new df
############
## get z-scores of just independent vars before going into the model
## You can get the zscore by using scale(x, center=TRUE, scale=TRUE)
## but need to make new df that has "parameter" which is the silicate or cover and "zscores" which is the standardized score
#percentcover_sef <- tibble("percentcover_zscore" = scale(SettlementDF_Full_C$percent_cover, center = TRUE, scale = TRUE))
#SettlementDF_Full_C$z_percent_cover <- scale(SettlementDF_Full_C$percent_cover, center = TRUE, scale = TRUE)
#SettlementDF_Full_C$z_SGD <- scale(SettlementDF_Full_C$silicate_avg, center = TRUE, scale = TRUE)
#SettlementDF_Full_C$z_SGD2 <- scale(SettlementDF_Full_C$silicate_avg^2, center = TRUE, scale = TRUE)
#SettlementDF_Full_V$z_percent_cover <- scale(SettlementDF_Full_V$percent_cover, center = TRUE, scale = TRUE)
#SettlementDF_Full_V$z_SGD <- scale(SettlementDF_Full_V$silicate_avg, center = TRUE, scale = TRUE)
#SettlementDF_Full_V$z_SGD2 <- scale(SettlementDF_Full_V$silicate_avg^2, center = TRUE, scale = TRUE)
### model the proportion of settlers that are alive out of the total on the plates
proportion_alive <- lm(proportion_alive ~ log(sum_total+1)*Location, data=SettlementDF_Full)
anova(proportion_alive) ## barely significant, but still significant, Location is not sig (as expected)
## Analysis of Variance Table
##
## Response: proportion_alive
## Df Sum Sq Mean Sq F value Pr(>F)
## log(sum_total + 1) 1 0.28120 0.281202 4.4035 0.04294 *
## Location 1 0.12179 0.121793 1.9072 0.17578
## log(sum_total + 1):Location 1 0.08614 0.086137 1.3489 0.25312
## Residuals 36 2.29892 0.063859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(proportion_alive) ## report R squared too in results
##
## Call:
## lm(formula = proportion_alive ~ log(sum_total + 1) * Location,
## data = SettlementDF_Full)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.56230 -0.14862 0.02246 0.13512 0.59074
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.76053 0.13328 5.706 1.72e-06 ***
## log(sum_total + 1) -0.10187 0.04996 -2.039 0.0489 *
## LocationVarari -0.32753 0.20191 -1.622 0.1135
## log(sum_total + 1):LocationVarari 0.08026 0.06911 1.161 0.2531
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2527 on 36 degrees of freedom
## Multiple R-squared: 0.1754, Adjusted R-squared: 0.1067
## F-statistic: 2.553 on 3 and 36 DF, p-value: 0.07069
check_model(proportion_alive)
### remove the 0s from both Varari and Cabral because skewing the data - 0s meaning the plates that had settlers, but none survived -- there are only a handful of recruits on those plates, so probably weren't even supposed to be there
### for all of the plates with 0, 6 or less settlers
proportion_alive_noZero <- lm(proportion_alive ~ log(sum_total+1)*Location, data=SettlementDF_Full %>%
filter(proportion_alive>0))
anova(proportion_alive_noZero) ## super significant, Location is not sig (as expected)
## Analysis of Variance Table
##
## Response: proportion_alive
## Df Sum Sq Mean Sq F value Pr(>F)
## log(sum_total + 1) 1 0.84218 0.84218 25.4635 1.738e-05 ***
## Location 1 0.00003 0.00003 0.0009 0.9765
## log(sum_total + 1):Location 1 0.01971 0.01971 0.5959 0.4458
## Residuals 32 1.05837 0.03307
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(proportion_alive_noZero)
##
## Call:
## lm(formula = proportion_alive ~ log(sum_total + 1) * Location,
## data = SettlementDF_Full %>% filter(proportion_alive > 0))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.309706 -0.089731 -0.006428 0.095317 0.307061
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.81691 0.09753 8.376 1.43e-09 ***
## log(sum_total + 1) -0.11284 0.03612 -3.124 0.00378 **
## LocationVarari 0.12431 0.17584 0.707 0.48473
## log(sum_total + 1):LocationVarari -0.04335 0.05616 -0.772 0.44582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1819 on 32 degrees of freedom
## Multiple R-squared: 0.4488, Adjusted R-squared: 0.3972
## F-statistic: 8.687 on 3 and 32 DF, p-value: 0.0002322
check_model(proportion_alive_noZero) ## high collinearity **
######################
## plot
######################
### with 0s
### don't worry about faceting by location because we are keeping both sites together in the model
densitydep_plot <- SettlementDF_Full %>%
ggplot(aes(x=log(sum_total+1),
y=proportion_alive)) +
geom_point(size=4, aes(color=Location)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4"),
limits=c("Cabral", "Varari")) +
geom_smooth(method="lm", color="black") +
theme_classic() +
labs(x= "Logged Total Settlers",
y= "Proportion of Settlers Alive") +
theme(axis.text.x=element_text(size=20),
axis.text.y=element_text(size=20),
axis.title.x=element_text(size=20),
axis.title.y=element_text(size=20),
legend.text=element_text(size=20), # Adjust legend text size
legend.title=element_text(size=20), # Adjust legend title size
strip.text=element_text(size=20))
densitydep_plot
## `geom_smooth()` using formula = 'y ~ x'
ggsave(here::here("Outputs", "RecruitmentTiles", "Thesis_Figures", "DensityPlot.jpg"),
width=12, height=9)
## `geom_smooth()` using formula = 'y ~ x'
### without 0s
densitydep_plot_noZeros <- SettlementDF_Full %>%
filter(proportion_alive>0) %>%
ggplot(aes(x=log(sum_total+1),
y=proportion_alive)) +
geom_point(size=4, aes(color=Location)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4"),
limits=c("Cabral", "Varari")) +
geom_smooth(method="lm", color="black") +
theme_classic() +
labs(x= "Logged Total Settlers",
y= "Proportion of Settlers Alive") +
theme(axis.text.x=element_text(size=20),
axis.text.y=element_text(size=20),
axis.title.x=element_text(size=20),
axis.title.y=element_text(size=20),
legend.text=element_text(size=20), # Adjust legend text size
legend.title=element_text(size=20), # Adjust legend title size
strip.text=element_text(size=20))
densitydep_plot_noZeros
## `geom_smooth()` using formula = 'y ~ x'
ggsave(here::here("Outputs", "RecruitmentTiles", "Thesis_Figures", "DensityPlot_NoZero.jpg"),
width=10, height=9)
## `geom_smooth()` using formula = 'y ~ x'
###########################
## student ts
###########################
############
## trying from stack overflow, old heavy package
#install_url('http://cran.r-project.org/src/contrib/Archive/heavy/heavy_0.38.19.tar.gz') ## not working right now
#fit <- heavyLm(m.marietta ~ CRSP, data = SettlementDF_Full, family = Student(df = 6))
#summary(fit)
#############
## trying with gaussian and identity link
#studentt_version1 <- stan_glm(proportion_alive ~ log(sum_total+1)*Location, data = SettlementDF_Full, family = gaussian(link = "identity"))
#summary(studentt_version1) ## not sure how to interpret
###########################
## quasibinomial
###########################
#RESULT: neither of these with "logit" link are significant
#BUT significant with "log" family. Without zero model did not converge.
#### with zeroes
quasibin_model <- glm(proportion_alive ~ log(sum_total+1)*Location, family = quasibinomial(link="log"), data=SettlementDF_Full)
summary(quasibin_model)
##
## Call:
## glm(formula = proportion_alive ~ log(sum_total + 1) * Location,
## family = quasibinomial(link = "log"), data = SettlementDF_Full)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.007783 0.198677 0.039 0.96897
## log(sum_total + 1) -0.294217 0.103082 -2.854 0.00711 **
## LocationVarari -0.851114 0.437535 -1.945 0.05959 .
## log(sum_total + 1):LocationVarari 0.241658 0.163563 1.477 0.14825
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.2676346)
##
## Null deviance: 14.483 on 39 degrees of freedom
## Residual deviance: 12.023 on 36 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 8
check_model(quasibin_model)
#### without zeroes
quasibin_model_noZero <- glm(proportion_alive ~ log(sum_total+1)*Location, family = quasibinomial(link="log"), data=SettlementDF_Full %>%
filter(proportion_alive>0))
## Warning: glm.fit: algorithm did not converge
summary(quasibin_model)
##
## Call:
## glm(formula = proportion_alive ~ log(sum_total + 1) * Location,
## family = quasibinomial(link = "log"), data = SettlementDF_Full)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.007783 0.198677 0.039 0.96897
## log(sum_total + 1) -0.294217 0.103082 -2.854 0.00711 **
## LocationVarari -0.851114 0.437535 -1.945 0.05959 .
## log(sum_total + 1):LocationVarari 0.241658 0.163563 1.477 0.14825
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.2676346)
##
## Null deviance: 14.483 on 39 degrees of freedom
## Residual deviance: 12.023 on 36 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 8
check_model(quasibin_model_noZero) ## some high collinearity
################################
## calculate residuals: plotting the residuals from previous model from hyp 2a against silicate
################################
#### with zeroes
resid(proportion_alive)
## 1 2 3 4 5 6
## 0.103180661 -0.234672191 -0.148615533 0.136257424 -0.284418902 -0.030158310
## 7 8 9 10 11 12
## 0.015331247 0.183737048 -0.067864920 -0.148615533 -0.211427814 0.282074671
## 13 14 15 16 17 18
## 0.310081414 0.351384467 0.351384467 -0.562304807 0.039411193 -0.161704446
## 19 20 21 22 23 24
## 0.029583750 0.047356116 -0.098284074 -0.004778701 -0.086284870 -0.235905358
## 25 26 27 28 29 30
## 0.114471662 0.590736950 -0.144692120 0.134742288 0.183355617 -0.145213635
## 31 32 33 34 35 36
## 0.205711886 0.124168843 -0.398227038 -0.409263050 0.009362440 0.058333295
## 37 38 39 40
## 0.080041717 0.378211732 -0.409263050 0.052775466
resid_proportionalive <- SettlementDF_Full %>%
bind_cols(residuals = resid(proportion_alive))
#### without zeroes
resid(proportion_alive_noZero)
## 1 2 3 4 5 6
## 0.084126854 -0.263781634 -0.192938565 0.128888055 -0.309705928 -0.057576059
## 7 8 9 10 11 12
## -0.011329355 0.153672717 -0.082967955 -0.192938565 -0.237380152 0.280461583
## 13 14 15 16 17 18
## 0.261308693 0.307061435 0.307061435 0.015344789 -0.193970993 -0.001526542
## 19 20 21 22 23 24
## 0.006190189 0.004963123 0.044612696 -0.045721926 -0.086965371 -0.098028525
## 25 26 27 28 29 30
## 0.230379314 -0.178309681 0.048519731 -0.044996476 -0.298249473 -0.061358086
## 31 32 33 34 35 36
## 0.231644210 -0.050385734 0.062440357 0.042406376 0.251306507 -0.052257043
resid_proportionalive_noZero <- SettlementDF_Full %>%
filter(proportion_alive>0) %>%
bind_cols(residuals = resid(proportion_alive_noZero))
#####################
## models
#####################
#### with zeroes
resids_silicateEffect_model <- lm(residuals ~ silicate_avg*Location, data=resid_proportionalive)
anova(resids_silicateEffect_model) ## very insignificant
## Analysis of Variance Table
##
## Response: residuals
## Df Sum Sq Mean Sq F value Pr(>F)
## silicate_avg 1 0.00012 0.000117 0.0019 0.9655
## Location 1 0.00008 0.000078 0.0013 0.9720
## silicate_avg:Location 1 0.07067 0.070671 1.1419 0.2924
## Residuals 36 2.22805 0.061890
summary(resids_silicateEffect_model)
##
## Call:
## lm(formula = residuals ~ silicate_avg * Location, data = resid_proportionalive)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53619 -0.17503 0.02615 0.13891 0.52805
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.09099 0.17254 0.527 0.601
## silicate_avg -0.02570 0.04612 -0.557 0.581
## LocationVarari -0.24427 0.24701 -0.989 0.329
## silicate_avg:LocationVarari 0.10373 0.09708 1.069 0.292
##
## Residual standard error: 0.2488 on 36 degrees of freedom
## Multiple R-squared: 0.03083, Adjusted R-squared: -0.04994
## F-statistic: 0.3817 on 3 and 36 DF, p-value: 0.7668
check_model(resids_silicateEffect_model)
#### without zeroes
resids_silicateEffect_model_noZero <- lm(residuals ~ silicate_avg*Location, data=resid_proportionalive_noZero)
anova(resids_silicateEffect_model_noZero)
## Analysis of Variance Table
##
## Response: residuals
## Df Sum Sq Mean Sq F value Pr(>F)
## silicate_avg 1 0.00093 0.000933 0.0286 0.8668
## Location 1 0.00053 0.000534 0.0163 0.8991
## silicate_avg:Location 1 0.01188 0.011880 0.3638 0.5507
## Residuals 32 1.04502 0.032657
summary(resids_silicateEffect_model_noZero)
##
## Call:
## lm(formula = residuals ~ silicate_avg * Location, data = resid_proportionalive_noZero)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33149 -0.07936 -0.01351 0.09740 0.31033
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.010988 0.126106 0.087 0.931
## silicate_avg -0.003151 0.034149 -0.092 0.927
## LocationVarari -0.092550 0.185599 -0.499 0.621
## silicate_avg:LocationVarari 0.044372 0.073567 0.603 0.551
##
## Residual standard error: 0.1807 on 32 degrees of freedom
## Multiple R-squared: 0.01261, Adjusted R-squared: -0.07996
## F-statistic: 0.1362 on 3 and 32 DF, p-value: 0.9377
check_model(resids_silicateEffect_model_noZero)
#####################
## plot
#####################
#### plot the residuals (this basically says that whatever points are above the best fit line are surviving more than expected due to density dependence because of SGD, anything below the line is surviving less than expected due to SGD)
#### with zeroes
survivorship_silicateEffect <- resid_proportionalive %>%
ggplot(aes(x=silicate_avg,
y=residuals)) +
geom_point(size=4, aes(color=Location)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4"),
limits=c("Cabral", "Varari")) +
geom_smooth(method="lm", color="black") +
theme_classic() +
labs(x = expression(Silicate~(mu*mol~L^-1)),
y= "Residuals of Alive Settlers") +
theme(axis.text.x=element_text(size=20),
axis.text.y=element_text(size=20),
axis.title.x=element_text(size=20),
axis.title.y=element_text(size=20),
legend.text=element_text(size=20), # Adjust legend text size
legend.title=element_text(size=20), # Adjust legend title size
strip.text=element_text(size=20))
survivorship_silicateEffect
## `geom_smooth()` using formula = 'y ~ x'
ggsave(here::here("Outputs", "RecruitmentTiles", "Thesis_Figures", "Survivorship_SilicateEffect.jpg"),
width=12, height=9)
## `geom_smooth()` using formula = 'y ~ x'
#### without zeroes
survivorship_silicateEffect_noZero <- resid_proportionalive_noZero %>%
ggplot(aes(x=silicate_avg,
y=residuals)) +
geom_point(size=4, aes(color=Location)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4"),
limits=c("Cabral", "Varari")) +
geom_smooth(method="lm", color="black") +
theme_classic() +
labs(x = expression(Silicate~(mu*mol~L^-1)),
y= "Residuals of Alive Settlers") +
theme(axis.text.x=element_text(size=20),
axis.text.y=element_text(size=20),
axis.title.x=element_text(size=20),
axis.title.y=element_text(size=20),
legend.text=element_text(size=20), # Adjust legend text size
legend.title=element_text(size=20), # Adjust legend title size
strip.text=element_text(size=20))
survivorship_silicateEffect_noZero
## `geom_smooth()` using formula = 'y ~ x'
ggsave(here::here("Outputs", "RecruitmentTiles", "Thesis_Figures", "Survivorship_SilicateEffect_NoZero.jpg"),
width=10, height=9)
## `geom_smooth()` using formula = 'y ~ x'
##################
## models
##################
## total alive
alive_predict_cover <- lm(percent_cover ~ sum_total_alive*Location, data =SettlementDF_Full)
anova(alive_predict_cover) ## insignificant, though logging makes it WAY better and almost significant (0.06)
## Analysis of Variance Table
##
## Response: percent_cover
## Df Sum Sq Mean Sq F value Pr(>F)
## sum_total_alive 1 70.53 70.525 2.2513 0.1422
## Location 1 76.39 76.390 2.4385 0.1271
## sum_total_alive:Location 1 83.55 83.554 2.6672 0.1112
## Residuals 36 1127.77 31.327
summary(alive_predict_cover)
##
## Call:
## lm(formula = percent_cover ~ sum_total_alive * Location, data = SettlementDF_Full)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.032 -3.120 -1.425 3.139 20.506
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.85949 1.43875 1.987 0.0545 .
## sum_total_alive 0.03253 0.07097 0.458 0.6495
## LocationVarari 0.36096 2.30200 0.157 0.8763
## sum_total_alive:LocationVarari 0.23165 0.14184 1.633 0.1112
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.597 on 36 degrees of freedom
## Multiple R-squared: 0.1697, Adjusted R-squared: 0.1005
## F-statistic: 2.452 on 3 and 36 DF, p-value: 0.07909
## proportion alive
propalive_predict_cover <- lm(percent_cover ~ proportion_alive*Location, data = SettlementDF_Full)
anova(propalive_predict_cover) ## significant
## Analysis of Variance Table
##
## Response: percent_cover
## Df Sum Sq Mean Sq F value Pr(>F)
## proportion_alive 1 142.31 142.314 4.6655 0.03751 *
## Location 1 33.99 33.994 1.1144 0.29816
## proportion_alive:Location 1 83.81 83.809 2.7475 0.10609
## Residuals 36 1098.12 30.503
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(propalive_predict_cover)
##
## Call:
## lm(formula = percent_cover ~ proportion_alive * Location, data = SettlementDF_Full)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.437 -3.247 -1.031 1.812 17.562
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.5028 2.7569 1.271 0.2120
## proportion_alive -0.6183 4.7914 -0.129 0.8980
## LocationVarari 6.9336 3.5295 1.964 0.0572 .
## proportion_alive:LocationVarari -11.4101 6.8836 -1.658 0.1061
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.523 on 36 degrees of freedom
## Multiple R-squared: 0.1915, Adjusted R-squared: 0.1241
## F-statistic: 2.843 on 3 and 36 DF, p-value: 0.05132
##################
## models without ZERO
##################
## total alive
alive_predict_cover_noZero <- lm(percent_cover ~ sum_total_alive*Location, data =SettlementDF_Full %>%
filter(percent_cover>0))
anova(alive_predict_cover_noZero) ## significant, better without logging
## Analysis of Variance Table
##
## Response: percent_cover
## Df Sum Sq Mean Sq F value Pr(>F)
## sum_total_alive 1 29.56 29.564 0.9619 0.33611
## Location 1 193.71 193.712 6.3024 0.01889 *
## sum_total_alive:Location 1 46.96 46.958 1.5278 0.22794
## Residuals 25 768.41 30.736
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(alive_predict_cover_noZero)
##
## Call:
## lm(formula = percent_cover ~ sum_total_alive * Location, data = SettlementDF_Full %>%
## filter(percent_cover > 0))
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8974 -2.7918 -0.9453 2.4895 17.7228
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.77664 1.61164 2.343 0.0274 *
## sum_total_alive 0.01777 0.07151 0.248 0.8058
## LocationVarari 2.88878 2.78900 1.036 0.3102
## sum_total_alive:LocationVarari 0.19126 0.15474 1.236 0.2279
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.544 on 25 degrees of freedom
## Multiple R-squared: 0.2602, Adjusted R-squared: 0.1714
## F-statistic: 2.931 on 3 and 25 DF, p-value: 0.05315
check_model(alive_predict_cover_noZero)
## proportion alive
propalive_predict_cover_noZero <- lm(percent_cover ~ proportion_alive*Location, data = SettlementDF_Full %>%
filter(percent_cover>0))
anova(propalive_predict_cover_noZero) ## very significant
## Analysis of Variance Table
##
## Response: percent_cover
## Df Sum Sq Mean Sq F value Pr(>F)
## proportion_alive 1 233.02 233.017 8.4472 0.007552 **
## Location 1 64.17 64.172 2.3263 0.139753
## proportion_alive:Location 1 51.83 51.827 1.8788 0.182651
## Residuals 25 689.63 27.585
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(propalive_predict_cover_noZero)
##
## Call:
## lm(formula = percent_cover ~ proportion_alive * Location, data = SettlementDF_Full %>%
## filter(percent_cover > 0))
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7592 -2.5642 -0.9592 2.0417 14.1464
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.917 3.267 1.811 0.0821 .
## proportion_alive -3.499 5.405 -0.647 0.5233
## LocationVarari 8.249 4.195 1.967 0.0604 .
## proportion_alive:LocationVarari -11.964 8.728 -1.371 0.1827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.252 on 25 degrees of freedom
## Multiple R-squared: 0.336, Adjusted R-squared: 0.2564
## F-statistic: 4.217 on 3 and 25 DF, p-value: 0.01523
check_model(propalive_predict_cover_noZero)
##################
## plots
##################
## total alive
alive_affect_cover <- SettlementDF_Full %>%
filter(percent_cover>0) %>% ## decide whether to keep in or not **
ggplot(aes(x=sum_total_alive,
y=percent_cover)) +
geom_point(size=4, aes(color=Location)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4"),
limits=c("Cabral", "Varari")) +
geom_smooth(method="lm", color="black") +
theme_classic() +
labs(x = "Total Alive Settlers",
y= "Percent Cover of P. acuta adults") +
theme(axis.text.x=element_text(size=20),
axis.text.y=element_text(size=20),
axis.title.x=element_text(size=20),
axis.title.y=element_text(size=20),
legend.text=element_text(size=20), # Adjust legend text size
legend.title=element_text(size=20), # Adjust legend title size
strip.text=element_text(size=20))
alive_affect_cover
## `geom_smooth()` using formula = 'y ~ x'
ggsave(here::here("Outputs", "RecruitmentTiles", "Thesis_Figures", "Cover_AliveEffect.jpg"),
width=12, height=9)
## `geom_smooth()` using formula = 'y ~ x'
## proportion alive
propalive_affect_cover <- SettlementDF_Full %>%
filter(percent_cover>0) %>%
ggplot(aes(x=proportion_alive,
y=percent_cover)) +
geom_point(size=4, aes(color=Location)) +
scale_color_manual(values=c("darkgoldenrod2", "firebrick4"),
limits=c("Cabral", "Varari")) +
geom_smooth(method="lm", color="black") +
theme_classic() +
labs(x = "Proportion of Alive Settlers",
y= "Percent Cover of P. acuta adults") +
theme(axis.text.x=element_text(size=20),
axis.text.y=element_text(size=20),
axis.title.x=element_text(size=20),
axis.title.y=element_text(size=20),
legend.text=element_text(size=20), # Adjust legend text size
legend.title=element_text(size=20), # Adjust legend title size
strip.text=element_text(size=20))
propalive_affect_cover
## `geom_smooth()` using formula = 'y ~ x'
ggsave(here::here("Outputs", "RecruitmentTiles", "Thesis_Figures", "Cover_PropAliveEffect.jpg"),
width=12, height=9)
## `geom_smooth()` using formula = 'y ~ x'